NETLOGO_PATH <- '/home/rick2/Apps/netlogo-5.1.0'
MODEL_PATH <- '/home/rick2/Documentos/engmat/hiv/codigo/netlogo/coutinho2.nlogo'
NSim <- 50
NTicks <- floor(12*(365/7))  # 625 ticks

tic <- Sys.time()

library(parallel)
processors <- detectCores()
cl <- makeCluster(processors)


# initialization
presim <- function(dummy, gui, nl.path, model.path) {
	library(RNetLogo)
	NLStart(nl.path, gui=gui)
	NLLoadModel(model.path)
}

# the simulation
simfun <- function(x, NTicks) {
	NLCommand('setup')
	d <- data.frame(T=numeric(NTicks), A=numeric(NTicks), D=numeric(NTicks))
	for(t in 1:NTicks) {
		T <- NLReport('count patches with [pcolor = T]')
		A <- NLReport('count patches with [pcolor = A1 or pcolor = A2]')
		D <- NLReport('count patches with [pcolor = D]')
		d[t,] <- c(T,A,D)
		NLCommand('update')
	}
	d
}

# quit
postsim <- function(x)
	NLQuit()


parLapply(cl, 1:processors, presim, gui=FALSE, nl.path=NETLOGO_PATH, model.path=MODEL_PATH)
result <- t(parSapply(cl, 1:NSim, simfun, NTicks))
parLapply(cl, 1:processors, postsim)
stopCluster(cl)

toc <- Sys.time()
print(toc-tic)

## Process results

Mean <- apply(result, 2, function(r) apply(matrix(unlist(r), nrow=NTicks), 1, mean))
Sdev <- apply(result, 2, function(r) apply(matrix(unlist(r), nrow=NTicks), 1, sd))
#MIN <- apply(result, 2, function(r) apply(matrix(unlist(r), nrow=NTicks), 1, min))
#MAX <- apply(result, 2, function(r) apply(matrix(unlist(r), nrow=NTicks), 1, max))

# standard error = standard deviation / sqrt(NSim)


library(reshape2)
Mean.short <- melt(data.frame(t=1:NTicks, Mean), id.vars='t')
Sdev.short <- melt(data.frame(t=1:NTicks, Sdev), id.vars='t')
Serr.short <- Sdev.short
Serr.short$value <- Serr.short$value / sqrt(NSim)
X <- data.frame(Mean.short, ymin=(Mean.short$value-Serr.short$value), ymax=(Mean.short$value+Serr.short$value))

# change time scales so that the first weeks are displayed with more promeninence
#Z$t[Z$t > 12] <- Z$t[Z$t > 12] + 100
#Z$t[Z$t <= 12] <- (Z$t[Z$t <= 12]-1) * 100/11

#output <- function(filename, width)
#	pdf(paste(filename, '.pdf', sep=''),
#		paper='special', width=8, height=8/width)

library(ggplot2)
#library(gridExtra)

theme_set(theme_classic() +
		theme(text=element_text('times')) +
		theme(strip.background=element_blank()) +
		theme(strip.text.x=element_text(face='bold',size=12)))

#output('virus', 2)
p <- ggplot(X, aes(x=t, y=value, color=variable, group=variable, ymin=ymin, ymax=ymax)) +
	geom_errorbar(colour='black', width=0.1) +
	geom_line() +
	geom_point(size=2, shape=21, fill='white') +
	scale_color_manual("Cell Type", values=c('green','orange','gray'))


print(p)

#Time difference of 8.284639 hours


	#geom_ribbon(aes(fill=variable), alpha=0.2) +
	#geom_line(aes(color=variable)) +
	#xlab("") + ylab(expression(paste('# cells x ', 10^5, sep=''))) +
	#scale_color_manual("Cell Type", values=c('green','yellow','red')) +
	#scale_fill_manual("Cell Type", values=c('green','yellow','red')) +
	#scale_x_continuous(breaks=c(1,50,100,seq(100+365/7,750,by=365/7)), labels=c(1,6,12,seq(2,13,by=1))) +
	#scale_y_continuous(breaks=seq(0,5e5,by=1e5), labels=seq(0,5,by=1)) +
	#geom_vline(xintercept=c(100), linetype="dotted") +
	#geom_text(data=data.frame(label=c("weeks","years"), x=c(50,400)), aes(label=label, x=x, y=y, group=NULL, ymin=NULL, ymax=NULL), y=-Inf, vjust=3.2, size=4) +
	#ggtitle("Cell count dynamics - no treatment")
#gt <- ggplot_gtable(ggplot_build(p))
#gt$layout$clip[gt$layout$name=="panel"] <- "off"
#grid.draw(gt)
#dev.off()

